home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / FOURN.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  3KB  |  86 lines
  1. PROCEDURE fourn(VAR data: gldarray; nn: glnnarray; ndim,isign: integer);
  2. (* Programs using routine FOURN must define the types
  3. TYPE
  4.    gldarray = ARRAY [1..ndat2] OF real;
  5.    glnnarray = ARRAY [1..ndim] OF integer;
  6. where ndat2 is twice the product of the transform lengths nn(i). *)
  7. VAR
  8.    i1,i2,i3,i2rev,i3rev,ibit,idim: integer;
  9.    ip1,ip2,ip3,ifp1,ifp2,k1,k2,n: integer;
  10.    ii1,ii2,ii3: integer;
  11.    nprev,nrem,ntot: integer;
  12.    tempi,tempr: real;
  13.    theta,wi,wpi,wpr,wr,wtemp: double;
  14. BEGIN
  15.    ntot := 1;
  16.    FOR idim := 1 TO ndim DO BEGIN
  17.       ntot := ntot*nn[idim]
  18.    END;
  19.    nprev := 1;
  20.    FOR idim := 1 TO ndim DO BEGIN
  21.       n := nn[idim];
  22.       nrem := ntot DIV (n*nprev);
  23.       ip1 := 2*nprev;
  24.       ip2 := ip1*n;
  25.       ip3 := ip2*nrem;
  26.       i2rev := 1;
  27.       FOR ii2 := 0 TO ((ip2-1) DIV ip1) DO BEGIN
  28.          i2 := 1+ii2*ip1;
  29.          IF (i2 < i2rev) THEN BEGIN
  30.             FOR ii1 := 0 TO ((ip1-2) DIV 2) DO BEGIN
  31.                i1 := i2+ii1*2;
  32.                FOR ii3 := 0 TO ((ip3-i1) DIV ip2) DO BEGIN
  33.                   i3 := i1+ii3*ip2;
  34.                   i3rev := i2rev+i3-i2;
  35.                   tempr := data[i3];
  36.                   tempi := data[i3+1];
  37.                   data[i3] := data[i3rev];
  38.                   data[i3+1] := data[i3rev+1];
  39.                   data[i3rev] := tempr;
  40.                   data[i3rev+1] := tempi
  41.                END
  42.             END
  43.          END;
  44.          ibit := ip2 DIV 2;
  45.          WHILE ((ibit >= ip1) AND (i2rev > ibit)) DO BEGIN
  46.             i2rev := i2rev-ibit;
  47.             ibit := ibit DIV 2
  48.          END;
  49.          i2rev := i2rev+ibit
  50.       END;
  51.       ifp1 := ip1;
  52.       WHILE (ifp1 < ip2) DO BEGIN
  53.          ifp2 := 2*ifp1;
  54.          theta := isign*6.28318530717959/(ifp2 DIV ip1);
  55.          wpr := -2.0*sqr(sin(0.5*theta));
  56.          wpi := sin(theta);
  57.          wr := 1.0;
  58.          wi := 0.0;
  59.          FOR ii3 := 0 TO ((ifp1-1) DIV ip1) DO BEGIN
  60.             i3 := 1+ii3*ip1;
  61.             FOR ii1 := 0 TO ((ip1-2) DIV 2) DO BEGIN
  62.                i1 := i3+ii1*2;
  63.                FOR ii2 := 0 TO ((ip3-i1) DIV ifp2) DO BEGIN
  64.                   i2 := i1+ii2*ifp2;
  65.                   k1 := i2;
  66.                   k2 := k1+ifp1;
  67.                   tempr := sngl(wr)*data[k2]
  68.                      -sngl(wi)*data[k2+1];
  69.                   tempi := sngl(wr)*data[k2+1]
  70.                      +sngl(wi)*data[k2];
  71.                   data[k2] := data[k1]-tempr;
  72.                   data[k2+1] := data[k1+1]-tempi;
  73.                   data[k1] := data[k1]+tempr;
  74.                   data[k1+1] := data[k1+1]+tempi
  75.                END
  76.             END;
  77.             wtemp := wr;
  78.             wr := wr*wpr-wi*wpi+wr;
  79.             wi := wi*wpr+wtemp*wpi+wi
  80.          END;
  81.          ifp1 := ifp2
  82.       END;
  83.       nprev := n*nprev
  84.    END
  85. END;
  86.